function distance3=distance3(x,data_par,Tradables_share,seed)
rng(seed);
Weight=[1/abs(data_par(1)),0,0;0,1/abs(data_par(2)),0;0,0,1/abs(data_par(3))];


%Vector of Initial Parameters
Par=[x(1),x(2),x(3)];


%%%%%First Set of Parameters%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%  Technical parameters  %%%%%%%%%%%%%%%%%%%%%%%%%
uptd=0.2; %Change to 0.1 if policy function fluctuate
outfreq=20; %Display frequency
iter_tol=600; % Iteration tolerance
tol=1e-5; %Numerical tolerance
 %Number of States
NB=100; %number of nodes.
bmin= -1;  %-1.11;
bmax=-0.;
NT      = 20;                % Number of grid points for tradables
NR      = 1;                % Number of grid points for nontradables
%NKAPPA  = 2;
NKAPPA  = 17;
NSS=NR*NT*NKAPPA;

%Moments from the Data
Debt_Data=data_par(1);
ProbabilityOFSS=data_par(2);
Std_CA=data_par(3);
ShareofT=Tradables_share;

%CALIBRATED PARAMETERS
beta=Par(1);
meank=Par(2);
sigmak=Par(3);

%rhok=Par(4);

%PARAMETERS NOT CALIBRATED
sigma=2;
ita=(1/0.83)-1;
rbar=0.04;
yN=1;

%AR(1) process for income
rho=0.460247378300448 ;
sigmay=0.032352010489055; 

rhok=rho;
 
%[Z,Zprob] = tauchenhussey(NT,0,rho,sigmay,sigmay); % Row are today, column are tomorrow
[Z,Zprob] = tauchen(NT,0,rho,sigmay,3); % Row are today, column are tomorrow

%[NewK,NewKprob]= tauchenhussey(NKAPPA,meank,rhok,sigmak,sigmak);
[NewK,NewKprob]= tauchen(NKAPPA,meank,rhok,sigmak,6);

Prob=kron(NewKprob,Zprob);   


yT      = exp(Z);
Sr     =  ones(NB,NR)*exp(rbar);
 

YT=repmat(yT,NKAPPA,NB )';
SR=repmat(Sr,1,NSS );
KAPPAS=kron(NewK,ones(NT,NB))';

%Calibration of OMEGA
%omega=0.3005; %Calibrated Outside of The Model
RatioOmega=mean((yT.*(1-ShareofT)).*(ShareofT.*(yT-Debt_Data*(exp(rbar)-1)).^(1+ita)).^(-1));
omega=1/(1+RatioOmega);

 %Repeated Stuff
 
 
B=zeros(NB,1);
B(1)=bmin;
 
for i=2:NB
    B(i)=B(i-1)+(bmax-B(i-1))/((NB-i+1)^1.05); 
end
 
 


%Q=1./SR;
q=1/exp(rbar);

b=repmat(B,1,NSS);
yn=ones(NB,NSS);
 
bp=  b;
%c= b./SR  +YT  -bp;
c= b  +YT  -bp*q;
%price=real((1-omega)/omega*c.^(1+ita));
price=real((1-omega)/omega*c.^(1+ita));

%c=ones(NB,NSS);
EE=zeros(NB,NSS);
cbind=c;
%bpmax=-kappa*(price+YT);
bpmax=-(1/q)*KAPPAS.*(price+YT);


%%%%%%%%%%%%%%%%%%DECENTRALIZED EQ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load A_de

iter=0;
d2=100;
%options = optimoptions('fsolve','Display','off');
disp('DE Iter      Norm');

  tol=0.00001
   iter=1
   
while d2>tol && iter <iter_tol    
    
oldp=price;
oldc=c;
oldbp=bp;

c_p=interp1(B,c,bp,'linear','extrap'); %This is 3D Array with consumptions tomorrow in today's grid



parfor i=1:NB
    for j=1:NSS
         c_tom=squeeze(c_p(i,j,:))'; % This is a vector with consumption tomorrow in each of the possible 15 states when today state is (i,j)
         totalcp=(omega*c_tom.^(-ita)+(1-omega)*yn(1,:).^(-ita));
         mu_p=real(omega*totalcp.^(sigma/ita).*(totalcp.^(-1/ita-1)).*(c_tom.^(-ita-1)));
         emu=beta*exp(rbar)*mu_p*Prob(j,:)'; %Expected Marginal Utility at state (i,j)
         EE(i,j)=(omega*cbind(i,j)^(-ita)+1-omega)^(sigma/ita -1/ita -1 )*omega*cbind(i,j)^(-ita-1)- emu;
%          EE(i,j)=(omega*cbind(i,j)^(-ita)+1-omega)^(sigma/ita -1/ita -1 )*omega*cbind(i,j)^(-ita-1)- emu(i,j);
        if EE(i,j)>1e-5
            bp(i,j)=bpmax(i,j);
            c(i,j)=cbind(i,j);
       
        else                            
            %f = @(cc) (omega*cc^(-ita)+1-omega)^(sigma/ita-1/ita-1)*omega*cc^(-ita-1)-emu;
            f = @(cc) ((1-omega)*cc^(ita)+omega)-(emu/omega)^(1/(sigma/ita-1/ita-1))*cc^((ita+1)/(sigma/ita-1/ita-1)+ita);
             lb=0;
             up=c(i,j)+1;
             x0=[lb, up];
             try
             [c(i,j), EE(i,j)]=fzero(f,x0);
             catch 
             x1=[lb, 2.5];
             [c(i,j), EE(i,j)]=fzero(f,x1);
             end                
%             x0=c(i,j);
%             [c(i,j), EE(i,j)]=fzero(f,x0);           
            bp(i,j)=max((1/q)*(YT(i,j)+b(i,j)-c(i,j)),bmin);
            bp(i,j)=min(bp(i,j),bmax);
        end
           
 
    end
end
 
price=real((1-omega)/omega*c.^(1+ita));

c=b+YT-max(q*bp,-KAPPAS.*(price+YT));  %qb'>=-kappa*(yt+yn*pn)
price=real((1-omega)/omega*c.^(1+ita));
bp=(1/q)*(b+YT-c);

d2=max([max(max(abs(c-oldc))),max(max(abs(bp-oldbp)))]);

%=====================Updating rule.Must be slow, important!==============
bp=uptd*bp+(1-uptd)*oldbp;
c=uptd*c+(1-uptd)*oldc;
price=uptd*price+(1-uptd)*oldp;
%========================================================================

bpmax=-(1/q)*KAPPAS.*(price.*yn+YT);
bpmax(bpmax>bmax)=bmax;
bpmax(bpmax<bmin)=bmin;
cbind=b+YT-q*bpmax;

 
cbind(cbind<0)=nan;


iter=iter+1;


D2(i,:)=d2;

metric(i)=d2;
if mod(iter, outfreq) == 0;
    fprintf('%d          %1.7f \n',iter,d2);
end
end
save de
clear  KAPPAS beta meank sigmak NewK NewKprob Prob Par data_par iter  
save A_de
clear all
load de

%%%%%%%%%%%%%%%%%%%%%%%%%%Simulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ISOM_01_model
%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%
%clear all
T=500000;
cut=100;  % Initial condition dependence
%load model_Brazil_k
%addpath('auxiliar')
clear kappa

%% 1. Simulation

[S_index,]=markov(Prob,T+cut+1,1,1:length(Prob));

    ySIM=YT(1,S_index)';  
    RSIM=SR(1,S_index)';  
    KAPPASIM=KAPPAS(1,S_index)';

bpSIM=zeros(T+cut,1);
%WfSIM=zeros(T+cut,1);
bindSIM=zeros(T+cut,1);
cSIM=zeros(T+cut,1);
pSIM=zeros(T+cut,1);

%bpSP_SIM=zeros(T+cut,1);
%bindSP_SIM=zeros(T+cut,1);
%cSP_SIM=zeros(T+cut,1);
%pSP_SIM=zeros(T+cut,1);

%Tau_SIM=zeros(T+cut,1);
%Tregion_SIM=zeros(T+cut,1);
 

bpSIM(1)=bmax+(bmin-bmax)/2;
%bpSP_SIM(1)=bmax+(bmin-bmax)/2;

for i=2:T+cut
     % de sim
    bpSIM(i)=interp1(B,bp(:,S_index(i)),bpSIM(i-1),'linear');
   cSIM(i)=bpSIM(i-1)+ySIM(i)-(1/(RSIM(i)))*bpSIM(i); 
 if isnan(bpSIM(i))==1
     %pause
 end
   % WfSIM(i)=interp1(B,Wf(:,S_index(i)),bpSIM(i-1),'linear');

    
         % sp sim
    %bpSP_SIM(i)=interp1(B,bpSP(:,S_index(i)),bpSP_SIM(i-1),'linear');
    %cSP_SIM(i)=bpSP_SIM(i-1)+ySIM(i)-(1/(RSIM(i)))*bpSP_SIM(i); 
    %Tau_SIM(i)=interp1(B,tau_sp(:,S_index(i)),bpSP_SIM(i-1),'linear');
end

pSIM=(1-omega)/omega*cSIM.^(1+ita); 
%pSP_SIM=(1-omega)/omega*cSP_SIM.^(1+ita);

bplim_SIM=-(RSIM).*  KAPPASIM.*(pSIM+ySIM);
%bplimSP_SIM=-(RSIM).*  KAPPASIM.*(pSP_SIM+ySIM);

%%%%% IMPORTANT TO USE TOLERANCE NOT <= %%%%%%%%%%%%
bindSIM((bpSIM-bplim_SIM)<=1e-5)=1;
%bindSP_SIM((bpSP_SIM-bplimSP_SIM)<=1e-5)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bpSIM(bindSIM==1)=bplim_SIM(bindSIM==1);
%bpSP_SIM(bindSP_SIM==1)=bplimSP_SIM(bindSP_SIM==1);

%Tregion_SIM(Tau_SIM>1)=1;

CA_SIM=q.*bpSIM(2:end)-bpSIM(1:end-1);  % >0 capital outflow
%CASP_SIM=bpSP_SIM(2:end)-bpSP_SIM(1:end-1);

% Eliminate Initial Condition Dependency
%S_index=S_index(cut+1:T+cut)';

bpSIM=bpSIM(cut+1:T+cut);
%bindSIM=bindSIM(cut+1:T+cut);
%bplim_SIM=bplim_SIM(cut+1:T+cut);
%cSIM=cSIM(cut+1:T+cut);
pSIM=pSIM(cut+1:T+cut);
CA_SIM=CA_SIM(cut:T+cut-1);

%bpSP_SIM=bpSP_SIM(cut+1:T+cut);
%bindSP_SIM=bindSP_SIM(cut+1:T+cut);
%bplimSP_SIM=bplimSP_SIM(cut+1:T+cut);
%cSP_SIM=cSP_SIM(cut+1:T+cut);
%pSP_SIM=pSP_SIM(cut+1:T+cut);
%CASP_SIM=CASP_SIM(cut:T+cut-1);

%WfSIM=WfSIM(cut+1:T+cut);
%Tau_SIM=Tau_SIM(cut+1:T+cut);
%Tregion_SIM=Tregion_SIM(cut+1:T+cut);
ySIM=ySIM(cut+1:T+cut);
%KAPPASIM=KAPPASIM(cut+1:T+cut);
%RSIM=RSIM(cut+1:T+cut);

Y_SIM=pSIM+ySIM;   % Income
%YSP_SIM=pSP_SIM+ySIM;

CA_SIM=CA_SIM./Y_SIM;  % Current Account
%CASP_SIM=CASP_SIM./YSP_SIM;

%CAchg_SIM=CA_SIM(2:end)-CA_SIM(1:end-1);
%CAchg_SIM=[0;CAchg_SIM];

Lev_SIM=bpSIM./Y_SIM;  % Leverage
%LevSP_SIM=bpSP_SIM./YSP_SIM;

%RER_SIM=(omega^(1/(1+ita))+(1-omega)^(1/(1+ita)).*pSIM.^(ita/(1+ita))).^((-1-ita)/ita); % Real Exchange Rate
%RERSP_SIM=(omega^(1/(1+ita))+(1-omega)^(1/(1+ita)).*pSP_SIM.^(ita/(1+ita))).^((-1-ita)/ita);


%MOMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear SS
SS=zeros(T,1);
SS(CA_SIM>(mean(CA_SIM)+2*std(CA_SIM)) )=1; % Crisis is defined as current account goes 2 sd away and collateral constraint binds in decentralized economy
MM1=mean(Lev_SIM);
MM2=mean(SS);
MM3=std(CA_SIM);
moments=[MM1,MM2,MM3];

% ca1=CA_SIM(1:end-1);
% ca2=CA_SIM(2:end);
% x_model=ones(2,length(ca1));
% x_model(2,:)=ca1;
% [breg,bregin,rreg,rintreg,stats]=regress(ca2,x_model');
% MM4=breg(2);
% MM3=stats(4);
%moments=[MM1,MM2,MM3,MM4];

distance3 =(data_par-moments)*Weight*(data_par-moments)';
end